#------------------------------------------------------------------------------
# Import libraries
#------------------------------------------------------------------------------

rm(list=ls())
library(LalRUtils)
libreq(tidyverse, data.table, zoo, tictoc, 
       fst, fixest, PanelMatch, patchwork,
       rio, magrittr, janitor, did, panelView, 
       ggiplot, tictoc, binsreg, interflex)
set.seed(42)
theme_set(lal_plot_theme())

#------------------------------------------------------------------------------





#------------------------------------------------------------------------------
# Define Paths
#------------------------------------------------------------------------------

# R studio
setwd( dirname(rstudioapi::getActiveDocumentContext()$path) )
# R default : unccoment if you use default R
# setwd(getSrcDirectory(function(){})[1])

#------------------------------------------------------------------------------




#------------------------------------------------------------------------------
# Load Data
#------------------------------------------------------------------------------

vcf <- fread("vcf_data_complete.csv", sep = ",")
setnames( vcf, "d", "D")
vcf_data <- copy( vcf )
gfc <- fread("gfc_dta.csv" )


#------------------------------------------------------------------------------



#------------------------------------------------------------------------------
# Figure 8: Pre-PESA Deforestation Rates by distance to mines (Panel A)
#------------------------------------------------------------------------------

vcf = vcf_data[year>= 1990]
vcf[, pref_bin := ntile(cover_1990, 10)]
vcf[, time_fra := year - 2008]
vcf[, D_fra := sch * (year >= 2008)]

# Generate subsample of data
xcreg = vcf[year == first_pesa_exposure - 1][min_dist_to_mine <= 1]
xcreg[, del_forest := cover_1990 - forest_index]
xcreg[, dist2 := min_dist_to_mine * 100]

# %% mining summary figure - binned scatterplot - fig8 panel A ;
gen = binsreg(y = xcreg$del_forest, 
              x = xcreg$dist2,
              nbins = 30, 
              polyreg = 2)
f = gen$bins_plot + lal_plot_theme() +
  labs(y = 'Decrease in Forest Index \nRelative to 1990 Baseline', 
       x = 'Distance (km)') +
  theme( legend.position = 'none', 
         axis.title = element_text( size = 22 ), 
         axis.text = element_text(size = 20))

# %%
ggsave("main_figure8PanelA.pdf", 
       f, width = 9, height = 5,
       device = cairo_pdf)


#------------------------------------------------------------------------------




#------------------------------------------------------------------------------
# Figure 8: Treatment effects on Annual Deforestation by Distance to Mines 
#           (Panel B)
#------------------------------------------------------------------------------


vcf_data <- copy( vcf )
vcf_data[, `:=`(id_f = as.factor(cellid),
                sty_f = as.factor(styear),
                year_f = as.factor(year),
                D_f    = as.factor(D)
)]
vcf_data[, distance := min_dist_to_mine * 110]
intsamp = vcf_data[distance <= 100]

#%% binning estimator - figure 8 panel B
mining_all = interflex(Y = "forest_index", 
                       D = "D_f", 
                       X = "distance",
                       FE = c("id_f", "sty_f"), 
                       cl = "id_f",
                       data = intsamp, 
                       theme.bw = T,
                       Xlabel = 'Distance', 
                       Dlabel = "PESA",
                       Ylabel = '\nForest Index (VCF)\n',
                       cutoffs = seq(0, 100, 10),
                       estimator = 'binning', CI = TRUE, 
                       cex.lab = 1.2,
                       cex.axis = 1.2, 
)

# %%
mining_all$figure 

ggsave("main_figure8PanelB.pdf", 
       mining_all$figure,
       width = 9, height = 5,
       device = cairo_pdf)


#------------------------------------------------------------------------------